from ‘Environmental enrichment alters protein expression as well as the proteomic response to cocaine in rat nucleus accumbens’ in [frontiers in behavioral neuroscience] https://doi.org/10.3389/fnbeh.2014.00246
Also see Zhang et al., 2016 from Green Group https://doi.org/10.1016/j.neuroscience.2016.09.051
Contents
library(readxl) # read_excel()
library(tidyverse) # dplyr, stringr, ggplot2, readr, tidyr
library(gt) # for gt tables
library(UniProt.ws) # 1st test - uniprot api
library(org.Rn.eg.db) # Rat genome annotation package - incl. AnnotationDbi
library(clusterProfiler) #gsea
library(enrichplot) # for dotplot()
# library(pathview) # commented out atm | (but for visual kegg mapping)
# library(msigdbr) # commented out atm | MSigDB gene sets
# # note some packages called directly instead e.g. cowplot:: & ggrepel::
df_1 <- read_excel("input/Lichti_2014_datasheet_1.XLSX",
sheet = "quant-all",
trim_ws = TRUE) # does not trim white space behind accessions?
colnames(df_1)
## [1] "Accession" "Name"
## [3] "Gene Symbol" "EC log2 fold"
## [5] "EC est fold" "EC p value"
## [7] "Coc log2 fold" "Coc est fold"
## [9] "Coc p value" "interaction log2 effect size"
## [11] "interaction effect size" "interaction p value"
## [13] "fraction"
str(df_1)
## tibble [2,835 × 13] (S3: tbl_df/tbl/data.frame)
## $ Accession : chr [1:2835] "Q3UHJ0" "Q3UHJ0" "P0C1X8" "P50475" ...
## $ Name : chr [1:2835] "AAK1_MOUSE" "AAK1_MOUSE" "AAK1_RAT" "SYAC_RAT" ...
## $ Gene Symbol : chr [1:2835] "AAK1" "AAK1" "Aak1" "AARS" ...
## $ EC log2 fold : num [1:2835] -0.442 0.404 0.216 -0.309 0.772 ...
## $ EC est fold : num [1:2835] -1.36 1.32 1.16 -1.24 1.71 ...
## $ EC p value : num [1:2835] 0.107 0.55 0.99 0.104 0.192 ...
## $ Coc log2 fold : num [1:2835] -0.312 -0.292 0.682 -0.15 -0.135 ...
## $ Coc est fold : num [1:2835] -1.24 -1.22 1.6 -1.11 -1.1 ...
## $ Coc p value : num [1:2835] 0.226 0.679 0.99 0.396 0.885 ...
## $ interaction log2 effect size: num [1:2835] 0.7 -0.357 0.195 0.442 -0.741 ...
## $ interaction effect size : num [1:2835] 1.62 -1.28 1.14 1.36 -1.67 ...
## $ interaction p value : num [1:2835] 0.0814 0.741 0.9899 0.1017 0.415 ...
## $ fraction : chr [1:2835] "cytosolic" "nuclear" "membrane" "cytosolic" ...
length(unique(df_1$Accession)) # 1915
## [1] 1915
head(df_1, 10)
## # A tibble: 10 × 13
## Accession Name `Gene Symbol` `EC log2 fold` `EC est fold` `EC p value`
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 Q3UHJ0 AAK1_MOUSE AAK1 -0.442 -1.36 0.107
## 2 Q3UHJ0 AAK1_MOUSE AAK1 0.404 1.32 0.550
## 3 P0C1X8 AAK1_RAT Aak1 0.216 1.16 0.990
## 4 P50475 SYAC_RAT AARS -0.309 -1.24 0.104
## 5 P50554 GABT_RAT ABAT 0.772 1.71 0.192
## 6 P50554 GABT_RAT ABAT 0.0918 1.07 0.389
## 7 P50554 GABT_RAT ABAT 0.670 1.59 0.990
## 8 P06795 MDR1B_MOUSE Abcb1b -2.35 -5.10 0.0906
## 9 Q5RKI8 ABCB8_RAT ABCB8 -0.875 -1.83 0.0632
## 10 Q5RKI8 ABCB8_RAT ABCB8 -0.0793 -1.06 0.991
## # ℹ 7 more variables: `Coc log2 fold` <dbl>, `Coc est fold` <dbl>,
## # `Coc p value` <dbl>, `interaction log2 effect size` <dbl>,
## # `interaction effect size` <dbl>, `interaction p value` <dbl>,
## # fraction <chr>
tail(df_1, 10)
## # A tibble: 10 × 13
## Accession Name `Gene Symbol` `EC log2 fold` `EC est fold` `EC p value`
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 P68511 1433F_RAT YWHAH -0.307 -1.24 0.0620
## 2 P68511 1433F_RAT YWHAH 0.263 1.20 0.605
## 3 P68511 1433F_RAT YWHAH 0.000514 1.00 1.000
## 4 P68254 1433T_MOUSE YWHAQ -0.220 -1.17 0.255
## 5 P68254 1433T_MOUSE YWHAQ -0.296 -1.23 0.947
## 6 P68254 1433T_MOUSE YWHAQ 0.491 1.41 0.990
## 7 P63101 1433Z_MOUSE YWHAZ -0.0697 -1.05 0.416
## 8 P63101 1433Z_MOUSE YWHAZ 0.470 1.39 0.479
## 9 P63101 1433Z_MOUSE YWHAZ 0.146 1.11 0.990
## 10 Q8VIL3 ZWINT_RAT Zwint -0.131 -1.09 0.541
## # ℹ 7 more variables: `Coc log2 fold` <dbl>, `Coc est fold` <dbl>,
## # `Coc p value` <dbl>, `interaction log2 effect size` <dbl>,
## # `interaction effect size` <dbl>, `interaction p value` <dbl>,
## # fraction <chr>
Getting info like protein names, etc
ratUp <- UniProt.ws(10116) # rat
my_ids <- df_1$Accession
# https://www.uniprot.org/help/return_fields
# https://www.uniprot.org/help/general_annotation
# takes the rat & mouse uniprot accession IDs in Lichti's df & maps further uniprot info onto them
up_rat_mouse <- mapUniProt(
from = "UniProtKB_AC-ID",
to = "UniProtKB",
columns = c("accession", "id", "reviewed", "protein_name", "gene_names",
"gene_primary", "organism_name", "length", "keywordid", "keyword", "xref_geneid"),
query = my_ids
#query = list(taxId = 10116, ids = my_ids) # isn't selective to only rats?
)
# add website links
# follows a certain stucture (doesn't account all but has most)
up_rat_mouse_entry <- up_rat_mouse %>%
mutate(uniprot_link = paste0("https://www.uniprot.org/uniprotkb/", Entry, "/entry")) %>% # add uniprot entry link!
mutate(wiki_link = paste0("https://en.wikipedia.org/wiki/", str_to_upper(Gene.Names..primary.))) # add wikipedia (needs 2 be uppercase)
#most of the wiki ones work e.g. in if 4 or more characters
# write out basic df
# write_excel_csv(up_rat_mouse_entry, "output/2025_12_08_lichti_2014_test_map_v1.csv")
df_2 <- dplyr::left_join(df_1, up_rat_mouse_entry,
by = join_by(Accession == From)) %>%
mutate(Direction = ifelse(`EC log2 fold` < 0, "Down", "Up"),
gene_name = str_to_title(`Gene.Names..primary.`),
.after = `Gene Symbol`) %>%
mutate(Significant = ifelse(`EC p value` <= 0.05, 'yes', 'no'), .after = `EC p value`) %>%
filter(str_detect(Name, "RAT")) # 2835 to 1506
## Warning in dplyr::left_join(df_1, up_rat_mouse_entry, by = join_by(Accession == : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 411 of `x` matches multiple rows in `y`.
## ℹ Row 1 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.
df_3 <- df_2 %>%
dplyr::select(Accession:Significant, fraction:wiki_link) %>%
mutate(across(where(is.numeric), ~ round(., 3))) %>% #if numeric col, round to 3 dp
group_by(fraction, Direction) %>%
arrange(`EC p value`, `EC log2 fold`)
# non significant plot - looks messy but you get the idea
ns_plot <- df_3 %>%
ggplot(aes(x = gene_name, y = `EC log2 fold`, color = Direction, alpha = Significant)) +
geom_point(shape = 1) +
geom_hline(yintercept = 0, color = 'grey20') +
scale_color_manual(breaks = c("Up", "Down"),
values = c("#dc322f", "#268bd2")) +
scale_alpha_manual(values = c(0.2,1)) +
labs(title = 'Lichti et al., 2014 - EC log2fc x gene name - NAc',
subtitle = 'Rats only (1506/2835), non-sig incl., split by fraction, LC-MS proteomics') +
facet_wrap(~fraction, space = 'free_x', scale = 'free_x') +
theme_bw()
ns_plot
# assign colors based on direction. | works by nvm ordering is weird
#label_colors <- ifelse(df_3$`EC log2 fold` > 0 , "#dc322f", "#268bd2")
# lollipop - ggplot ver - a big ugly
sig_plot <- df_3 %>%
filter(`EC p value` <= 0.05) %>% # 139
ggplot(aes(x = reorder(gene_name, `EC log2 fold`), # orders by log fold
y = `EC log2 fold`,
color = Direction)) +
geom_segment(aes(x = gene_name,
yend = `EC log2 fold`, y = 0)) +
geom_point(aes(alpha = `EC p value`),
shape = 19, size = 4) +
scale_alpha(range = c(1,.1), guide = 'none') + # rev alpha so darker = more significant
# geom_text(aes(label = gene_name), color = "grey10", size = 2.5,
# position = position_stack()) +
geom_hline(yintercept = 0, color = 'grey20', linewidth = 0.1) +
scale_color_manual(breaks = c("Up", "Down"),
values = c("#dc322f", "#268bd2"),
guide = 'none') +
labs(title = 'Lichti et al., 2014 - EC Log2FC x Gene Name - NAc',
subtitle = 'Rats only (139/1506/2835), sig only, by fraction, LC-MS proteomics',
y = 'EC Log2FC',
x = 'Gene Name') +
#coord_flip() + #fliptheplot
facet_wrap(~fraction, space = 'free_x', scale = 'free_x') +
#theme_void() +
theme_bw() +
theme(
#axis.title = element_blank(),
#axis.text.x = element_blank()
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, # 90, hj1 t/R, vj.5 cent
size = 9, face = 'italic'),
axis.text.y = element_text(hjust = 1, color = "grey20", size = 10), # rota 90 & hjust=1 (top/right)
)
sig_plot
volcano_adjpv <- ggplot(df_3,
aes(x =`EC log2 fold`,
y = -log2(`EC p value`), label = gene_name)) +
geom_point(aes(color = case_when(
`EC p value` < 0.05 & `EC log2 fold` > 0 ~ "Upregulated (p<0.05)",
`EC p value` < 0.05 & `EC log2 fold` < 0 ~ "Downregulated (p<0.05)",
`EC p value` < 0.1 & `EC log2 fold` > 0 ~ "Upregulated (p<0.1)",
`EC p value` < 0.1 & `EC log2 fold` < 0 ~ "Downregulated (p<0.1)",
TRUE ~ "Not Significant")),
alpha = 1, size = 2.5, shape = 1) +
xlim(-7.5,7.5) +
geom_vline(xintercept = -1.5, linewidth =0.2, linetype = "dashed", color = "grey50") +
geom_vline(xintercept = 1.5, linewidth =0.2, linetype = "dashed", color = "grey50") +
geom_hline(yintercept = -log2(0.05), linetype = "dashed", linewidth = 0.2, color = "grey30") +
scale_color_manual(values = c("firebrick3",
"royalblue4",
"lightpink1",
"lightblue3" ,
"grey80"),
breaks = c("Upregulated (p<0.05)",
"Downregulated (p<0.05)",
"Upregulated (p<0.1)",
"Downregulated (p<0.1)",
"Not Significant"),
labels = c("Upregulated Proteins (up to p<0.05)",
"Downregulated Proteins (up to p<0.05)",
"Upregulated Proteins (up to p<0.1)",
"Downregulated Proteins (up to p<0.1)" ,
"Not Significant")) +
ggrepel::geom_text_repel() +
facet_wrap(~ fraction, ncol = 1, scales = "free_y") +
labs(
title = "Lichti et al., 2014 LC-MS proteomics NAc rats - main effect of housing (EE vs IC) (2025_12_26)",
subtitle = "Split by fraction | Volcano plots made from their provided df, not 100% sure on their p-value if adj or not?",
x = "log2FC",
y = "-Log2 (p-value)",
color = "Protein Regulation") +
geom_hline(yintercept = -log2(0.05), linetype = "dashed", linewidth = 0.2, color = "grey30") +
# annotate("text", x = .86, y = -log2(0.05), label = "Adj. p-value cutoff (0.05) \n",
# fontface = 'italic', size = 4, color = "grey10") +
# geom_hline(yintercept = -log2(0.1), linetype = "dashed", linewidth = 0.2, color = "grey50") +
# annotate("text", x = .86, y = -log2(0.1), label = "Adj. p-value cutoff (0.1) \n",
# fontface = 'italic', size = 4, color = "grey30") +
theme_bw() +
theme(legend.position = "bottom",
strip.text = element_text(size = 12, face = "bold", color = "grey10")
)
volcano_adjpv
## Warning: ggrepel: 626 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 331 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 441 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
gt table for sig proteins only
df_3 |>
dplyr::select(fraction, gene_name, `Protein.names`, Direction,
`EC log2 fold`, `EC p value`, Entry) |> #cols for table
dplyr::filter(`EC p value` <= 0.05) |>
gt(
groupname_col = dplyr::group_vars(df_3)
) |>
tab_header(
title = "Lichti et al., 2014 - EC vs. IC - NAc",
subtitle = "Cleaned lc-ms proteomics dataframe, signficant rats only (139/1506/2835
), organised by fraction + direction & arranged by pvalue + EC log2FC"
) |>
opt_stylize(style = 1, color = 'gray', add_row_striping = TRUE)
| Lichti et al., 2014 - EC vs. IC - NAc | ||||
| Cleaned lc-ms proteomics dataframe, signficant rats only (139/1506/2835 ), organised by fraction + direction & arranged by pvalue + EC log2FC | ||||
| gene_name | Protein.names | EC log2 fold | EC p value | Entry |
|---|---|---|---|---|
| cytosolic - Up | ||||
| Glod4 | Glyoxalase domain-containing protein 4 | 0.538 | 0.007 | Q5I0D1 |
| Gsto1 | Glutathione S-transferase omega-1 (GSTO-1) (EC 2.5.1.18) (Glutathione S-transferase omega 1-1) (GSTO 1-1) (Glutathione-dependent dehydroascorbate reductase) (EC 1.8.5.1) (Monomethylarsonic acid reductase) (MMA(V) reductase) (EC 1.20.4.2) (S-(Phenacyl)glutathione reductase) (SPG-R) | 0.389 | 0.013 | Q9Z339 |
| Tpi1 | Triosephosphate isomerase (TIM) (EC 5.3.1.1) (Methylglyoxal synthase) (EC 4.2.3.3) (Triose-phosphate isomerase) | 0.274 | 0.016 | P48500 |
| Prdx5 | Peroxiredoxin-5, mitochondrial (EC 1.11.1.24) (Antioxidant enzyme B166) (AOEB166) (PLP) (Peroxiredoxin V) (Prx-V) (Peroxisomal antioxidant enzyme) (Thioredoxin peroxidase PMP20) (Thioredoxin-dependent peroxiredoxin 5) | 0.293 | 0.018 | Q9R063 |
| Arpc5l | Actin-related protein 2/3 complex subunit 5-like protein (Arp2/3 complex 16 kDa subunit 2) (ARC16-2) | 0.510 | 0.018 | A1L108 |
| Wdr7 | WD repeat-containing protein 7 (TGF-beta resistance-associated protein TRAG) | 0.652 | 0.018 | Q9ERH3 |
| Akr1a1 | Aldo-keto reductase family 1 member A1 (EC 1.1.1.2) (EC 1.1.1.372) (EC 1.1.1.54) (3-DG-reducing enzyme) (Alcohol dehydrogenase [NADP(+)]) (Aldehyde reductase) (Glucuronate reductase) (EC 1.1.1.19) (Glucuronolactone reductase) (EC 1.1.1.20) (S-nitroso-CoA reductase) (ScorR) (EC 1.6.-.-) | 0.258 | 0.020 | P51635 |
| Fahd2a | Oxaloacetate tautomerase Fahd2a, mitochondrial (EC 5.3.2.2) (Fumarylacetoacetate hydrolase domain-containing protein 2A) | 0.365 | 0.020 | B2RYW9 |
| Slc25a11 | Mitochondrial 2-oxoglutarate/malate carrier protein (OGCP) (alpha-oxoglutarate carrier) (Solute carrier family 25 member 11) (SLC25A11) | 0.430 | 0.020 | P97700 |
| Kyat3 | Kynurenine--oxoglutarate transaminase 3 (EC 2.6.1.7) (Cysteine-S-conjugate beta-lyase 2) (EC 4.4.1.13) (Kynurenine aminotransferase 3) (Kynurenine aminotransferase III) (KATIII) (Kynurenine--glyoxylate transaminase) (EC 2.6.1.63) (Kynurenine--oxoglutarate transaminase III) | 0.680 | 0.020 | Q58FK9 |
| Atp5f1b | ATP synthase F(1) complex catalytic subunit beta, mitochondrial (EC 7.1.2.2) (ATP synthase F1 subunit beta) | 0.161 | 0.021 | P10719 |
| Dnm1 | Dynamin-1 (EC 3.6.5.5) (B-dynamin) (D100) (Dynamin I) (Dynamin, brain) | 0.249 | 0.021 | P21575 |
| Dlat | Dihydrolipoyllysine-residue acetyltransferase component of pyruvate dehydrogenase complex, mitochondrial (EC 2.3.1.12) (70 kDa mitochondrial autoantigen of primary biliary cirrhosis) (PBC) (Dihydrolipoamide acetyltransferase component of pyruvate dehydrogenase complex) (Pyruvate dehydrogenase complex component E2) (PDC-E2) (PDCE2) | 0.273 | 0.021 | P08461 |
| Stip1 | Stress-induced-phosphoprotein 1 (STI1) (Hsc70/Hsp90-organizing protein) (Hop) | 0.299 | 0.021 | O35814 |
| Qdpr | Dihydropteridine reductase (EC 1.5.1.34) (HDHPR) (Quinoid dihydropteridine reductase) | 0.351 | 0.021 | P11348 |
| Septin3 | Neuronal-specific septin-3 (G-septin) (P40) | 0.571 | 0.021 | Q9WU34 |
| Slc4a10 | Sodium-driven chloride bicarbonate exchanger (Solute carrier family 4 member 10) | 3.018 | 0.021 | Q80ZA5 |
| Crmp1 | Dihydropyrimidinase-related protein 1 (DRP-1) (Collapsin response mediator protein 1) (CRMP-1) (Inactive dihydropyrimidinase) | 0.190 | 0.022 | Q62950 |
| Atp5po | ATP synthase peripheral stalk subunit OSCP, mitochondrial (ATP synthase subunit O) (Oligomycin sensitivity conferral protein) (OSCP) (Sperm flagella protein 4) | 0.394 | 0.022 | Q06647 |
| Nudt3 | Diphosphoinositol polyphosphate phosphohydrolase 1 (DIPP-1) (EC 3.6.1.52) (Diadenosine hexaphosphate hydrolase) (Ap6A hydrolase) (EC 3.6.1.61) (Endopolyphosphatase) (EC 3.6.1.10) (Nucleoside diphosphate-linked moiety X motif 3) (Nudix motif 3) (m7GpppN-mRNA hydrolase) (EC 3.6.1.62) (m7GpppX diphosphatase) (EC 3.6.1.59) | 0.793 | 0.025 | Q566C7 |
| Fahd1 | Oxaloacetate tautomerase FAHD1, mitochondrial (EC 5.3.2.2) (Acylpyruvase FAHD1) (EC 3.7.1.5) (Fumarylacetoacetate hydrolase domain-containing protein 1) (Oxaloacetate decarboxylase) (OAA decarboxylase) (EC 4.1.1.112) | 1.525 | 0.027 | Q6AYQ8 |
| Uba1 | Ubiquitin-like modifier-activating enzyme 1 (EC 6.2.1.45) (Ubiquitin-activating enzyme E1) | 0.197 | 0.028 | Q5U300 |
| Atp5f1c | ATP synthase F(1) complex subunit gamma, mitochondrial (ATP synthase F1 subunit gamma) (F-ATPase gamma subunit) | 0.384 | 0.028 | P35435 |
| Phgdh | D-3-phosphoglycerate dehydrogenase (3-PGDH) (EC 1.1.1.95) | 0.450 | 0.028 | O08651 |
| Mpst | 3-mercaptopyruvate sulfurtransferase (MST) (EC 2.8.1.2) | 0.332 | 0.029 | P97532 |
| Hibadh | 3-hydroxyisobutyrate dehydrogenase, mitochondrial (HIBADH) (EC 1.1.1.31) | 0.409 | 0.029 | P29266 |
| Nme1 | Nucleoside diphosphate kinase A (NDK A) (NDP kinase A) (EC 2.7.4.6) (Metastasis inhibition factor NM23) (Tumor metastatic process-associated protein) | 0.299 | 0.030 | Q05982 |
| Cryab | Alpha-crystallin B chain (Alpha(B)-crystallin) | 0.621 | 0.030 | P23928 |
| Hspa4 | Heat shock 70 kDa protein 4 (Ischemia responsive 94 kDa protein) | 0.200 | 0.031 | O88600 |
| Cs | Citrate synthase, mitochondrial (EC 2.3.3.1) (Citrate (Si)-synthase) | 0.322 | 0.031 | Q8VHF5 |
| Hspa8 | Heat shock cognate 71 kDa protein (EC 3.6.4.10) (Heat shock 70 kDa protein 8) | 0.208 | 0.032 | P63018 |
| Aco1 | Cytoplasmic aconitate hydratase (Aconitase) (EC 4.2.1.3) (Citrate hydro-lyase) (Iron regulatory protein 1) (IRP1) (Iron-responsive element-binding protein 1) (IRE-BP 1) | 0.271 | 0.032 | Q63270 |
| Atp5mf | ATP synthase F(0) complex subunit f, mitochondrial (ATP synthase membrane subunit f) | 0.332 | 0.032 | D3ZAF6 |
| Slc12a5 | Solute carrier family 12 member 5 (Electroneutral potassium-chloride cotransporter 2) (Furosemide-sensitive K-Cl cotransporter) (K-Cl cotransporter 2) (rKCC2) (Neuronal K-Cl cotransporter) | 0.370 | 0.032 | Q63633 |
| Slc6a17 | Sodium-dependent neutral amino acid transporter SLC6A17 (Sodium-dependent neurotransmitter transporter NTT4) (Solute carrier family 6 member 17) | 0.379 | 0.032 | P31662 |
| Dpp7 | Dipeptidyl peptidase 2 (EC 3.4.14.2) (Dipeptidyl aminopeptidase II) (Dipeptidyl peptidase 7) (Dipeptidyl peptidase II) (DPP II) (Quiescent cell proline dipeptidase) | 0.467 | 0.032 | Q9EPB1 |
| Atp5mg | ATP synthase F(0) complex subunit g, mitochondrial (ATPase subunit g) (ATP synthase membrane subunit g) | 0.599 | 0.032 | Q6PDU7 |
| Actrt1 | Actin-related protein T1 | 1.752 | 0.032 | Q5XIK1 |
| Hsp90aa1 | Heat shock protein HSP 90-alpha (EC 3.6.4.10) (Heat shock 86 kDa) (HSP 86) (HSP86) | 0.336 | 0.033 | P82995 |
| Coq9 | Ubiquinone biosynthesis protein COQ9, mitochondrial | 1.045 | 0.033 | Q68FT1 |
| Dlg2 | Disks large homolog 2 (Channel-associated protein of synapse-110) (Chapsyn-110) (Postsynaptic density protein PSD-93) | 1.625 | 0.033 | Q63622 |
| Exoc2 | Exocyst complex component 2 (Exocyst complex component Sec5) (rSec5) | 1.769 | 0.033 | O54921 |
| Ppp5c | Serine/threonine-protein phosphatase 5 (PP5) (EC 3.1.3.16) (Protein phosphatase T) (PPT) | 2.036 | 0.033 | P53042 |
| Mpc2 | Mitochondrial pyruvate carrier 2 (Brain protein 44) (Protein 0-44) | 0.736 | 0.035 | P38718 |
| Ctsd | Cathepsin D (EC 3.4.23.5) [Cleaved into: Cathepsin D 12 kDa light chain; Cathepsin D 9 kDa light chain; Cathepsin D 34 kDa heavy chain; Cathepsin D 30 kDa heavy chain] | 1.318 | 0.035 | P24268 |
| Nfasc | Neurofascin | 0.172 | 0.036 | P97685 |
| Gnao1 | Guanine nucleotide-binding protein G(o) subunit alpha (EC 3.6.5.-) | 0.196 | 0.036 | P59215 |
| Pygb | Glycogen phosphorylase, brain form (EC 2.4.1.1) | 0.229 | 0.036 | P53534 |
| Tubb2a | Tubulin beta-2A chain | 0.316 | 0.036 | P85108 |
| Ppp2r2a | Serine/threonine-protein phosphatase 2A 55 kDa regulatory subunit B alpha isoform (PP2A subunit B isoform B55-alpha) (B55) (PP2A subunit B isoform BRA) (PP2A subunit B isoform PR55-alpha) (PP2A subunit B isoform R2-alpha) (PP2A subunit B isoform alpha) | 0.321 | 0.036 | P36876 |
| Ran | GTP-binding nuclear protein Ran (EC 3.6.5.-) (GTPase Ran) (Ras-like protein TC4) (Ras-related nuclear protein) | 0.505 | 0.036 | P62828 |
| Dpp10 | Inactive dipeptidyl peptidase 10 (Dipeptidyl peptidase X) (DPP X) (Kv4 potassium channel auxiliary subunit) | 0.579 | 0.036 | Q6Q629 |
| Vcp | Transitional endoplasmic reticulum ATPase (TER ATPase) (EC 3.6.4.6) (15S Mg(2+)-ATPase p97 subunit) (Valosin-containing protein) (VCP) | 0.184 | 0.037 | P46462 |
| Idh3b | Isocitrate dehydrogenase [NAD] subunit beta, mitochondrial (Isocitric dehydrogenase subunit beta) (NAD(+)-specific ICDH subunit beta) | 0.236 | 0.037 | Q68FX0 |
| Tubb3 | Tubulin beta-3 chain (Neuron-specific class III beta-tubulin) | 0.616 | 0.040 | Q4QRB4 |
| Usp7 | Ubiquitin carboxyl-terminal hydrolase 7 (EC 3.4.19.12) (Deubiquitinating enzyme 7) (Herpesvirus-associated ubiquitin-specific protease) (rHAUSP) (Ubiquitin thioesterase 7) (Ubiquitin-specific-processing protease 7) | 0.914 | 0.040 | Q4VSI4 |
| Mdh1 | Malate dehydrogenase, cytoplasmic (EC 1.1.1.37) (Aromatic alpha-keto acid reductase) (KAR) (EC 1.1.1.96) (Cytosolic malate dehydrogenase) | 0.279 | 0.041 | O88989 |
| Pgrmc1 | Membrane-associated progesterone receptor component 1 (mPR) (25-DX) (Acidic 25 kDa protein) (Ventral midline antigen) (VEMA) | 0.281 | 0.041 | P70580 |
| Fah | Fumarylacetoacetase (FAA) (EC 3.7.1.2) (Beta-diketonase) (Fumarylacetoacetate hydrolase) | 0.303 | 0.041 | P25093 |
| Idh1 | Isocitrate dehydrogenase [NADP] cytoplasmic (IDH) (IDH1) (EC 1.1.1.42) (Cytosolic NADP-isocitrate dehydrogenase) (IDPc) (NADP(+)-specific ICDH) (Oxalosuccinate decarboxylase) | 0.331 | 0.041 | P41562 |
| Cstb | Cystatin-B (Cystatin-beta) (Liver thiol proteinase inhibitor) (Stefin-B) | 0.444 | 0.041 | P01041 |
| Capza1 | F-actin-capping protein subunit alpha-1 (CapZ alpha-1) | 0.568 | 0.041 | B2GUZ5 |
| Slc3a2 | Amino acid transporter heavy chain SLC3A2 (4F2 cell-surface antigen heavy chain) (4F2hc) (Solute carrier family 3 member 2) (CD antigen CD98) | 0.188 | 0.042 | Q794F9 |
| Sncb | Beta-synuclein (Phosphoneuroprotein 14) (PNP 14) | 0.315 | 0.042 | Q63754 |
| Ola1 | Obg-like ATPase 1 | 0.316 | 0.042 | A0JPJ7 |
| Lap3 | Cytosol aminopeptidase (EC 3.4.11.1) (Cysteinylglycine-S-conjugate dipeptidase) (EC 3.4.13.23) (Leucine aminopeptidase 3) (LAP-3) (Leucyl aminopeptidase) (LAP) (Peptidase S) (Proline aminopeptidase) (EC 3.4.11.5) (Prolyl aminopeptidase) | 0.474 | 0.042 | Q68FS4 |
| Inpp4a | Inositol polyphosphate-4-phosphatase type I A (Inositol polyphosphate 4-phosphatase type I) (Type I inositol 3,4-bisphosphate 4-phosphatase) (EC 3.1.3.66) | 0.225 | 0.043 | Q62784 |
| Pak1 | Serine/threonine-protein kinase PAK 1 (EC 2.7.11.1) (Alpha-PAK) (Protein kinase MUK2) (p21-activated kinase 1) (PAK-1) (p68-PAK) | 0.294 | 0.044 | P35465 |
| Gapdh | Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (EC 1.2.1.12) (38 kDa BFA-dependent ADP-ribosylation substrate) (BARS-38) (Peptidyl-cysteine S-nitrosylase GAPDH) (EC 2.6.99.-) | 0.378 | 0.044 | P04797 |
| Cmpk1 | UMP-CMP kinase (EC 2.7.4.14) (Deoxycytidylate kinase) (CK) (dCMP kinase) (Nucleoside-diphosphate kinase) (EC 2.7.4.6) (Uridine monophosphate/cytidine monophosphate kinase) (UMP/CMP kinase) (UMP/CMPK) | 0.462 | 0.044 | Q4KM73 |
| Tkt | Transketolase (TK) (EC 2.2.1.1) | 0.192 | 0.045 | P50137 |
| Gna11 | Guanine nucleotide-binding protein subunit alpha-11 (G alpha-11) (G-protein subunit alpha-11) | 0.275 | 0.045 | Q9JID2 |
| Gpi | Glucose-6-phosphate isomerase (GPI) (EC 5.3.1.9) (Autocrine motility factor) (AMF) (Neuroleukin) (NLK) (Phosphoglucose isomerase) (PGI) (Phosphohexose isomerase) (PHI) | 0.161 | 0.046 | Q6P6V0 |
| Got1 | Aspartate aminotransferase, cytoplasmic (cAspAT) (EC 2.6.1.1) (EC 2.6.1.3) (Cysteine aminotransferase, cytoplasmic) (Cysteine transaminase, cytoplasmic) (cCAT) (Glutamate oxaloacetate transaminase 1) (Transaminase A) | 0.290 | 0.046 | P13221 |
| Psmb2 | Proteasome subunit beta type-2 (Macropain subunit C7-I) (Multicatalytic endopeptidase complex subunit C7-I) (Proteasome component C7-I) (Proteasome subunit beta-4) (beta-4) | 0.695 | 0.046 | P40307 |
| Tcp1 | T-complex protein 1 subunit alpha (TCP-1-alpha) (EC 3.6.1.-) (CCT-alpha) | 0.202 | 0.048 | P28480 |
| Tuba1b | Tubulin alpha-1B chain (EC 3.6.5.-) (Alpha-tubulin 2) (Tubulin alpha-2 chain) [Cleaved into: Detyrosinated tubulin alpha-1B chain] | 2.398 | 0.048 | Q6P9V9 |
| Sod2 | Superoxide dismutase [Mn], mitochondrial (EC 1.15.1.1) | 0.914 | 0.049 | P07895 |
| Akr1b1 | Aldo-keto reductase family 1 member B1 (EC 1.1.1.21) (EC 1.1.1.300) (EC 1.1.1.372) (EC 1.1.1.54) (Aldehyde reductase) (Aldose reductase) (AR) | 0.189 | 0.050 | P07943 |
| Pacsin1 | Protein kinase C and casein kinase substrate in neurons protein 1 (Dynamin proline-rich domain-interacting protein) (Dynamin PRD-interacting protein) (Synaptic, dynamin-associated protein I) (Syndapin-1) (Syndapin-I) (SdpI) | 0.189 | 0.050 | Q9Z0W5 |
| Prdx6 | Peroxiredoxin-6 (EC 1.11.1.27) (1-Cys peroxiredoxin) (1-Cys PRX) (Acidic calcium-independent phospholipase A2) (aiPLA2) (EC 3.1.1.4) (Antioxidant protein 2) (Glutathione-dependent peroxiredoxin) (Lysophosphatidylcholine acyltransferase 5) (LPC acyltransferase 5) (LPCAT-5) (Lyso-PC acyltransferase 5) (EC 2.3.1.23) (Non-selenium glutathione peroxidase) (NSGPx) (Thiol-specific antioxidant protein) | 0.195 | 0.050 | O35244 |
| Sfxn3 | Sideroflexin-3 | 0.272 | 0.050 | Q9JHY2 |
| Bdh1 | D-beta-hydroxybutyrate dehydrogenase, mitochondrial (EC 1.1.1.30) (3-hydroxybutyrate dehydrogenase) (BDH) | 0.343 | 0.050 | P29147 |
| Hnrnpd | Heterogeneous nuclear ribonucleoprotein D0 (hnRNP D0) (AU-rich element RNA-binding protein 1) | 0.856 | 0.050 | Q9JJ54 |
| cytosolic - Down | ||||
| Arf4 | ADP-ribosylation factor 4 | -2.212 | 0.013 | P61751 |
| Adgrl1 | Adhesion G protein-coupled receptor L1 (Calcium-independent alpha-latrotoxin receptor 1) (CIRL-1) (Latrophilin-1) | -3.426 | 0.014 | O88917 |
| C1qbp | Complement component 1 Q subcomponent-binding protein, mitochondrial (GC1q-R protein) (Glycoprotein gC1qBP) (C1qBP) | -2.686 | 0.014 | O35796 |
| Cygb | Cytoglobin (Nitric oxygen dioxygenase CYGB) (EC 1.14.12.-) (Nitrite reductase CYGB) (EC 1.7.-.-) (Pseudoperoxidase CYGB) (EC 1.11.1.-) (Stellate cell activation-associated protein) (Superoxide dismutase CYGB) (EC 1.15.1.1) | -3.009 | 0.018 | Q921A4 |
| Gyg1 | Glycogenin-1 (GN-1) (GN1) (EC 2.4.1.186) | -2.965 | 0.018 | O08730 |
| Entpd2 | Ectonucleoside triphosphate diphosphohydrolase 2 (NTPDase 2) (EC 3.6.1.-) (CD39 antigen-like 1) (Ecto-ATP diphosphohydrolase 2) (Ecto-ATPDase 2) (Ecto-ATPase 2) | -2.420 | 0.018 | O35795 |
| Map2k2 | Dual specificity mitogen-activated protein kinase kinase 2 (MAP kinase kinase 2) (MAPKK 2) (EC 2.7.12.2) (ERK activator kinase 2) (MAPK/ERK kinase 2) (MEK 2) | -2.246 | 0.018 | P36506 |
| Vps35l | VPS35 endosomal protein-sorting factor-like | -0.910 | 0.018 | Q5XI83 |
| Csad | Cysteine sulfinic acid decarboxylase (EC 4.1.1.29) (Aspartate 1-decarboxylase) (EC 4.1.1.11) (Cysteine-sulfinate decarboxylase) (Sulfinoalanine decarboxylase) | -5.304 | 0.020 | Q64611 |
| Psmb7 | Proteasome subunit beta type-7 (EC 3.4.25.1) (Macropain chain Z) (Multicatalytic endopeptidase complex chain Z) (Proteasome subunit Z) (Proteasome subunit beta-2) (beta-2) | -3.307 | 0.020 | Q9JHW0 |
| Pnpo | Pyridoxine-5'-phosphate oxidase (EC 1.4.3.5) (Pyridoxamine-phosphate oxidase) | -3.138 | 0.020 | O88794 |
| Rabggta | Geranylgeranyl transferase type-2 subunit alpha (EC 2.5.1.60) (Geranylgeranyl transferase type II subunit alpha) (Rab geranyl-geranyltransferase subunit alpha) (Rab GG transferase alpha) (Rab GGTase alpha) (Rab geranylgeranyltransferase subunit alpha) | -1.157 | 0.020 | Q08602 |
| Ddx1 | ATP-dependent RNA helicase DDX1 (EC 3.6.4.13) (DEAD box protein 1) | -0.622 | 0.020 | Q641Y8 |
| Copb1 | Coatomer subunit beta (Beta-coat protein) (Beta-COP) | -0.455 | 0.020 | P23514 |
| Stx1a | Syntaxin-1A (Neuron-specific antigen HPC-1) (Synaptotagmin-associated 35 kDa protein) (P35A) | -0.431 | 0.025 | P32851 |
| Serpina1 | Alpha-1-antiproteinase (Alpha-1-antitrypsin) (Alpha-1-proteinase inhibitor) (Serpin A1) | -2.763 | 0.028 | P17475 |
| Fkbp1a | Peptidyl-prolyl cis-trans isomerase FKBP1A (PPIase FKBP1A) (EC 5.2.1.8) (12 kDa FK506-binding protein) (12 kDa FKBP) (FKBP-12) (FK506-binding protein 1A) (FKBP-1A) (Immunophilin FKBP12) (Rotamase) | -0.952 | 0.030 | Q62658 |
| Adk | Adenosine kinase (AK) (EC 2.7.1.20) (Adenosine 5'-phosphotransferase) | -2.509 | 0.031 | Q64640 |
| Man2c1 | Alpha-mannosidase 2C1 (EC 3.2.1.24) (Alpha-D-mannoside mannohydrolase) (Mannosidase alpha class 2C member 1) | -2.007 | 0.031 | P21139 |
| Por | NADPH--cytochrome P450 reductase (CPR) (P450R) (EC 1.6.2.4) | -1.797 | 0.031 | P00388 |
| Lcmt1 | Leucine carboxyl methyltransferase 1 (EC 2.1.1.233) ([Phosphatase 2A protein]-leucine-carboxy methyltransferase 1) | -1.726 | 0.031 | Q6P4Z6 |
| Ppm1e | Protein phosphatase 1E (EC 3.1.3.16) (Ca(2+)/calmodulin-dependent protein kinase phosphatase N) (CaMKP-N) (CaMKP-nucleus) (CaMKN) (Partner of PIX 1) (Partner of PIX-alpha) (Partner of PIXA) | -1.311 | 0.031 | Q80Z30 |
| Pafah1b2 | Platelet-activating factor acetylhydrolase IB subunit alpha2 (EC 3.1.1.47) (PAF acetylhydrolase 30 kDa subunit) (PAF-AH 30 kDa subunit) (PAF-AH subunit beta) (PAFAH subunit beta) (Platelet-activating factor acetylhydrolase alpha 2 subunit) (PAF-AH alpha 2) | -0.666 | 0.031 | O35264 |
| Rpn2 | Dolichyl-diphosphooligosaccharide--protein glycosyltransferase subunit 2 (Dolichyl-diphosphooligosaccharide--protein glycosyltransferase 63 kDa subunit) (Ribophorin II) (RPN-II) (Ribophorin-2) | -0.681 | 0.033 | P25235 |
| Rala | Ras-related protein Ral-A (EC 3.6.5.2) | -1.185 | 0.035 | P63322 |
| Prepl | Prolyl endopeptidase-like (EC 3.4.21.-) (Prolylendopeptidase-like) | -0.682 | 0.035 | Q5HZA6 |
| Plcb1 | 1-phosphatidylinositol 4,5-bisphosphate phosphodiesterase beta-1 (EC 3.1.4.11) (PLC-154) (Phosphoinositide phospholipase C-beta-1) (Phospholipase C-I) (PLC-I) (Phospholipase C-beta-1) (PLC-beta-1) | -0.253 | 0.039 | P10687 |
| Prune1 | Exopolyphosphatase PRUNE1 (EC 3.6.1.1) | -0.557 | 0.041 | Q6AYG3 |
| Slc8a1 | Sodium/calcium exchanger 1 (Na(+)/Ca(2+)-exchange protein 1) (Solute carrier family 8 member 1) | -0.490 | 0.041 | Q01728 |
| Nono | Non-POU domain-containing octamer-binding protein (NonO protein) | -2.858 | 0.043 | Q5FVM4 |
| Pir | Pirin (EC 1.13.11.24) (Probable quercetin 2,3-dioxygenase PIR) (Probable quercetinase) | -0.444 | 0.045 | Q5M827 |
| Sacm1l | Phosphatidylinositol-3-phosphatase SAC1 (EC 3.1.3.64) (Phosphatidylinositol-4-phosphate phosphatase) (Suppressor of actin mutations 1-like protein) | -1.268 | 0.050 | Q9ES21 |
| Ndufv2 | NADH dehydrogenase [ubiquinone] flavoprotein 2, mitochondrial (EC 7.1.1.2) (NADH-ubiquinone oxidoreductase 24 kDa subunit) | -0.608 | 0.050 | P19234 |
| nuclear - Up | ||||
| Atp5me | ATP synthase F(0) complex subunit e, mitochondrial (ATPase subunit e) (ATP synthase membrane subunit e) | 1.618 | 0.017 | P29419 |
| Phb1 | Prohibitin 1 | 1.098 | 0.022 | P67779 |
| Uqcrc1 | Cytochrome b-c1 complex subunit 1, mitochondrial (Complex III subunit 1) (Core protein I) (Ubiquinol-cytochrome-c reductase complex core protein 1) | 0.853 | 0.023 | Q68FY0 |
| Atp5po | ATP synthase peripheral stalk subunit OSCP, mitochondrial (ATP synthase subunit O) (Oligomycin sensitivity conferral protein) (OSCP) (Sperm flagella protein 4) | 0.953 | 0.023 | Q06647 |
| Timm13 | Mitochondrial import inner membrane translocase subunit Tim13 | 0.708 | 0.034 | P62076 |
| Ndufa9 | NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 9, mitochondrial (Complex I-39kD) (CI-39kD) (NADH-ubiquinone oxidoreductase 39 kDa subunit) (Sperm flagella protein 3) | 0.561 | 0.045 | Q5BK63 |
| Atp5f1a | ATP synthase F(1) complex subunit alpha, mitochondrial (ATP synthase F1 subunit alpha) | 0.748 | 0.045 | P15999 |
| Phb2 | Prohibitin-2 (B-cell receptor-associated protein BAP37) (BAP-37) | 0.751 | 0.045 | Q5XIH7 |
| Uqcrc2 | Cytochrome b-c1 complex subunit 2, mitochondrial (Complex III subunit 2) (Core protein II) (Ubiquinol-cytochrome-c reductase complex core protein 2) | 0.909 | 0.045 | P32551 |
| Ndufv2 | NADH dehydrogenase [ubiquinone] flavoprotein 2, mitochondrial (EC 7.1.1.2) (NADH-ubiquinone oxidoreductase 24 kDa subunit) | 0.948 | 0.045 | P19234 |
| Uqcrfs1 | Cytochrome b-c1 complex subunit Rieske, mitochondrial (EC 7.1.1.8) (Complex III subunit 5) (Cytochrome b-c1 complex subunit 5) (Liver regeneration-related protein LRRGT00195) (Rieske iron-sulfur protein) (RISP) (Rieske protein UQCRFS1) (Ubiquinol-cytochrome c reductase iron-sulfur subunit) [Cleaved into: Cytochrome b-c1 complex subunit 9 (Su9) (Subunit 9) (8 kDa subunit 9) (Complex III subunit IX) (Cytochrome b-c1 complex subunit 11) (UQCRFS1 mitochondrial targeting sequence) (UQCRFS1 MTS) (Ubiquinol-cytochrome c reductase 8 kDa protein)] | 0.989 | 0.045 | P20788 |
| Ppia | Peptidyl-prolyl cis-trans isomerase A (PPIase A) (EC 5.2.1.8) (Cyclophilin A) (Cyclosporin A-binding protein) (Rotamase A) (p1B15) (p31) [Cleaved into: Peptidyl-prolyl cis-trans isomerase A, N-terminally processed] | 1.167 | 0.045 | P10111 |
| Cd47 | Leukocyte surface antigen CD47 (Integrin-associated protein) (IAP) (CD antigen CD47) | 0.538 | 0.047 | P97829 |
| Stxbp1 | Syntaxin-binding protein 1 (N-Sec1) (Protein unc-18 homolog 1) (Unc18-1) (Protein unc-18 homolog A) (Unc-18A) (p67) (rbSec1) | 0.663 | 0.047 | P61765 |
| Aldoa | Fructose-bisphosphate aldolase A (EC 4.1.2.13) (Muscle-type aldolase) | 0.675 | 0.047 | P05065 |
| Tomm22 | Mitochondrial import receptor subunit TOM22 homolog (rTOM22) (Translocase of outer membrane 22 kDa subunit homolog) | 0.796 | 0.047 | Q75Q41 |
| Maoa | Amine oxidase [flavin-containing] A (EC 1.4.3.21) (EC 1.4.3.4) (Monoamine oxidase type A) (MAO-A) | 0.818 | 0.047 | P21396 |
| Cox5a | Cytochrome c oxidase subunit 5A, mitochondrial (Cytochrome c oxidase polypeptide Va) | 0.833 | 0.047 | P11240 |
| Cs | Citrate synthase, mitochondrial (EC 2.3.3.1) (Citrate (Si)-synthase) | 0.900 | 0.047 | Q8VHF5 |
| Pkm | Pyruvate kinase PKM (EC 2.7.1.40) (Pyruvate kinase muscle isozyme) (Threonine-protein kinase PKM2) (EC 2.7.11.1) (Tyrosine-protein kinase PKM2) (EC 2.7.10.2) | 0.908 | 0.047 | P11980 |
| Atp5f1b | ATP synthase F(1) complex catalytic subunit beta, mitochondrial (EC 7.1.2.2) (ATP synthase F1 subunit beta) | 0.930 | 0.047 | P10719 |
| Aldoc | Fructose-bisphosphate aldolase C (EC 4.1.2.13) (Brain-type aldolase) | 1.140 | 0.047 | P09117 |
| Hbb | Hemoglobin subunit beta-1 (Beta-1-globin) (Hemoglobin beta chain, major-form) (Hemoglobin beta-1 chain) | 1.286 | 0.047 | P02091 |
| Cisd1 | CDGSH iron-sulfur domain-containing protein 1 (Cysteine transaminase CISD1) (EC 2.6.1.3) (MitoNEET) | 1.344 | 0.047 | B0K020 |
| Arhgdia | Rho GDP-dissociation inhibitor 1 (Rho GDI 1) (Rho-GDI alpha) | 1.765 | 0.047 | Q5XI73 |
| Hsph1 | Heat shock protein 105 kDa (Heat shock 110 kDa protein) | 0.427 | 0.050 | Q66HA8 |
| nuclear - Down | ||||
| Dlg4 | Disks large homolog 4 (Postsynaptic density protein 95) (PSD-95) (Synapse-associated protein 90) (SAP-90) (SAP90) | -0.378 | 0.047 | P31016 |
| Shank3 | SH3 and multiple ankyrin repeat domains protein 3 (Shank3) (Proline-rich synapse-associated protein 2) (ProSAP2) (SPANK-2) | -0.439 | 0.049 | Q9JLU4 |
already partially prepped from above
# # ranked combined list
# geneList_df <- df_3 %>%
# mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
# group_by(Entry) %>%
# summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
# arrange(-ranking)
#
# # convert 2 a named vector to supply to fgsea()
# geneList <- tibble::deframe(geneList_df) # 955
#
# plot(geneList)
# # ranked cytosolic list
# geneList_df_c <- df_3 %>%
# filter(fraction == 'cytosolic') %>%
# mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
# group_by(Entry) %>%
# summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
# arrange(-ranking)
#
# # convert 2 a named vector to supply to fgsea()
# geneList_c <- tibble::deframe(geneList_df_c) # 689
#
# plot(geneList_c)
# # ranked membrane list
# geneList_df_m <- df_3 %>%
# filter(fraction == 'membrane') %>%
# mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
# group_by(Entry) %>%
# summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
# arrange(-ranking)
#
# # convert 2 a named vector to supply to fgsea()
# geneList_m <- tibble::deframe(geneList_df_m) # 340
#
# plot(geneList_m)
# # ranked nuclear list
# geneList_df_n <- df_3 %>%
# filter(fraction == 'nuclear') %>%
# mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
# group_by(Entry) %>%
# summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
# arrange(-ranking)
#
# # convert 2 a named vector to supply to fgsea()
# geneList_n <- tibble::deframe(geneList_df_n) # 477
#
# plot(geneList_n)
# # nuclear
# kegg_res_n <- gseKEGG(geneList = geneList_n,
# organism = 'rno',
# keyType = 'uniprot',
# by = 'fgsea',
# pvalueCutoff = 0.05,
# pAdjustMethod = 'BH', #BejaminiHoc Correciton
# )
#
# gse_kegg_tib_n <- as_tibble(kegg_res_n) # 21 sig
# enrichplot::dotplot(kegg_res_n,
# split = ".sign", # split by + vs - enrirchment score
# showCategory = 30,
# title = "Nuclear - NAc - Lichti et al., 2014 (EE vs IC) - gseKEGG_r - adj.p<0.05") +
# facet_grid(.~.sign) # split plot by seign)
# # cytosolic
# kegg_res_c <- gseKEGG(geneList = geneList_c,
# organism = 'rno',
# keyType = 'uniprot',
# by = 'fgsea',
# pvalueCutoff = 0.05,
# pAdjustMethod = 'BH', #BejaminiHoc Correciton
# )
#
# gse_kegg_tib_c <- as_tibble(kegg_res_c) # 21 sig
# enrichplot::dotplot(kegg_res_c,
# split = ".sign", # split by + vs - enrirchment score
# showCategory = 30,
# title = "Cytosolic - NAc - Lichti et al., 2014 (EE vs IC) - gseKEGG_r - adj.p<0.05") +
# facet_grid(.~.sign) # split plot by seign)
# # membrane
# kegg_res_m <- gseKEGG(geneList = geneList_m,
# organism = 'rno',
# keyType = 'uniprot',
# by = 'fgsea',
# pvalueCutoff = 0.05,
# pAdjustMethod = 'BH', #BejaminiHoc Correciton
# )
#
# gse_kegg_tib_m <- as_tibble(kegg_res_m) # 3 sig
# enrichplot::dotplot(kegg_res_m,
# split = ".sign", # split by + vs - enrirchment score
# showCategory = 30,
# title = "Membrane - NAc - Lichti et al., 2014 (EE vs IC) - gseKEGG_r - adj.p<0.05") +
# facet_grid(.~.sign) # split plot by seign)
# # oxidative phosphorylation - rno00190
# rno00190 <- pathview(gene.data = geneList_n,
# gene.idtype = "SYMBOL",
# pathway.id = 'rno00190', # oxidative phosphorylation - Rat
# species = 'rno',
# limit = list(gene=max(abs(geneList)), cpd=1))
# make vector of fractions (cytosolic, membrane, & nuclear) - in prep of sequences
fraction_names <- unique(df_3$fraction)
# make container lists to store results
# gsea_results_list <- vector(mode = "list", length = 10)
gsea_results_list <- vector(mode = "list") # actually make empty
for (fr in fraction_names){ # SEQUENCE - for each item (fr) in fraction_names
# OPERATIONS
# Part 1
# make a ranked gene list for this fraciton
geneList_working_df <- df_3 %>%
# filter by item | fraction is column in df_3, fr is i (current item in seq)
filter(fraction == fr) %>%
# compute ranking metric
mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
# for any duplicate entries, # average & remove any nas
group_by(Entry) %>%
summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
# lastly, arrange from pos to neg
arrange(-ranking)
# store in list (not necessary...)
gsea_results_list$geneList_df[[fr]] <- geneList_working_df
# Part 1.2
# next, deframe into a named vector needed to supply to fgsea()
geneList_vec_working <- tibble::deframe(geneList_working_df)
# store in list
gsea_results_list$geneList_vec[[fr]] <- geneList_vec_working
# Part 1.3 | plots fine but does not save?!?!?!
# Plot ranked gene list
gL_plot <- plot(geneList_vec_working, main = paste(fr, '- Ranked gene list'))
print(gL_plot)
gsea_results_list$geneList_plots[[fr]] <- gL_plot
# Part 2
# run KEGG gsea with clusterProfiler
kegg_res_working <- clusterProfiler::gseKEGG(geneList = geneList_vec_working,
organism = 'rno',
keyType = 'uniprot',
by = 'fgsea',
pvalueCutoff = 0.05,
pAdjustMethod = 'BH', #BejaminiHoc Correciton
)
# save kegg result original to list vec
gsea_results_list$kegg_results[[fr]] <- kegg_res_working
# not needed here so compute & save directly? or can do anyway? idk?
gsea_results_list$kegg_results_tib[[fr]] <- tibble::as_tibble(kegg_res_working)
# Part 3
# Make dotplots | this one plots & saves properly!
# 3.1 - regular dot blot
dp <- enrichplot::dotplot(kegg_res_working,
showCategory = 30,
font.size = 11,
title = paste0(stringr::str_glue(
'{fr} - NAc - Lichti, 2014 (EE vs IC) - gseKEGG_r')))
print(dp) # print dot plots
gsea_results_list$dotplots[[fr]] <- dp # save to list
# 3.2 - dot plot split by sign
dp_sign <- enrichplot::dotplot(kegg_res_working,
split = '.sign',
showCategory = 30,
font.size = 11,
title = paste0(stringr::str_glue(
'{fr} - NAc - Lichti, 2014 (EE vs IC) - gseKEGG_r split by sign'))) +
facet_grid(.~.sign)
print(dp_sign) # print dot plots
gsea_results_list$dotplots_sign[[fr]] <- dp_sign # save to list
}
## NULL
## Reading KEGG annotation online: "https://rest.kegg.jp/link/rno/pathway"...
## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/rno"...
## Reading KEGG annotation online: "https://rest.kegg.jp/conv/uniprot/rno"...
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## leading edge analysis...
## done...
## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...
## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## leading edge analysis...
## done...
# horrific name but at least clear...
#Lichti_2014_gsea_results_2025_12_12_v1 <- gsea_results_list
# save(Lichti_2014_gsea_results_2025_12_12_v1,
# file = "output/2025_12_12_Lichti_2014_gsea_KEGG_results_rats_NAc_v1.rda")
# Warning means that some genes have tied rankings
# note the membrane one is especially high 89.71% - this corresponds to the weird # pattern in geneList ranking plot -
# Warning messages:
# 1: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
# There are ties in the preranked stats (33.82% of the list).
# The order of those tied genes will be arbitrary, which may produce unexpected results.
# 2: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
# There are ties in the preranked stats (30.61% of the list).
# The order of those tied genes will be arbitrary, which may produce unexpected results.
# 3: In fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize, :
# For some of the pathways the P-values were likely overestimated. For such pathways log2err is set to NA.
# 4: In fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize, :
# For some pathways, in reality P-values are less than 1e-10. You can set the `eps` argument to zero for better estimation.
# 5: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
# There are ties in the preranked stats (89.71% of the list).
# The order of those tied genes will be arbitrary, which may produce unexpected results.
check_1 <- gsea_results_list$geneList_df$cytosolic
length(unique(gsea_results_list$geneList_df$cytosolic$Entry)) # 689 unique proteins
## [1] 689
length(unique(gsea_results_list$geneList_df$cytosolic$ranking)) # 456 unique ranks
## [1] 456
# cow plot for combined + plotlist = to plot a list of plots
cowplot::plot_grid(plotlist = c(gsea_results_list$dotplots,
gsea_results_list$dotplots_sign),
ncol = 2)
# print(gsea_results_list$dotplots)
# weird... highlights importance of plotting data!!!
plot(df_3$`EC log2 fold`)
2025_12_12 - Note to self maybe add some testing & progress of loop info - ref to epi cook book. Also try nested for loops with GO analysis e.g. for each of the 3 fractions (cytosolic, nuclear, membrane) run each of the 3 GO (CC, BP, MF)
colnames(df_3) # check col names
## [1] "Accession" "Name" "Gene Symbol"
## [4] "Direction" "gene_name" "EC log2 fold"
## [7] "EC est fold" "EC p value" "Significant"
## [10] "fraction" "Entry" "Entry.Name"
## [13] "Reviewed" "Protein.names" "Gene.Names"
## [16] "Gene.Names..primary." "Organism" "Length"
## [19] "Keyword.ID" "Keywords" "GeneID"
## [22] "uniprot_link" "wiki_link"
# lit review table
lit_re <- df_3 %>%
mutate(
Space = " ", # empty space column
Year = 2014,
Author = 'Lichti',
Animal = 'Rat',
Condition = paste0("Fraction - ", fraction),
Region = 'Nucleus_Accumbens',
Whole = 'No',
Exp_type = 'Protein',
Method = 'PROTEOMICS_LCMS_LFQ_OS',
Used_name = `Gene Symbol`,
Mod = " ",
Mod_type = " ",
Gene_name = gene_name, # lazy way to avoid moving rows (don't need headers anyway)
Comparison = 1,
In_EE = Direction, # note was coded as >0 Up, <0 Down. Maybe change one to same??
Fold_change = `EC log2 fold`,
#Arrows = ifelse(Direction == 'Up', '+', '-'),
Arrows = case_when( # dplyr case when, instead of ifelse
`EC log2 fold` > 2.5 ~ '+++',
`EC log2 fold` > 1 ~ '++',
`EC log2 fold` > .05 ~ '+',
`EC log2 fold` > -.05 ~ '=',
`EC log2 fold` > -1 ~ '-',
`EC log2 fold` > -2.5 ~ '--',
`EC log2 fold` <= -2.5 ~ '---',
),
#Significant = ifelse(`EC p value` <= 0.05, 'yes', 'no'),
signficant = Significant,
pvalue = `EC p value`,
N_tested = " ",
Stats_adj = 'adjusted',
Note = ' ',
Full_name = Protein.names,
Uniprot_rat_accession = Entry,
Uniprot_entry_name = Entry.Name,
Uniprot_link = uniprot_link,
Wikipedia_link = wiki_link,
Curation_method = ifelse(!is.na(Entry), 'auto - uniprot_acc_query', 'manual'),
Class_group = ' ',
Other_link = ' ',
Other_comments = '*misc col shows the col names in their excel df',
misc = paste0(Accession, ";", `Gene Symbol`, ";", Name)
)
#write out | na values = blank
# write_excel_csv(lit_re,
# "output/2025_12_12_Lichti_2014_rats_NAc_data_lit_review_v1.csv",
# na = "")
2025_12_24 - while most of the reanalysis matches the results presented in the paper. some did not quite. I suspected this was due to using the mouse database (which i excluded from my analysis). The code below is checking that.
# list of EE vs IC sig proteins shown in table 1 of the paper | 15 genes/proteins
acc_2_check <- c('P61751', 'Q64611', 'P24268', 'Q641Y8', 'Q63622', 'Q6PDL0',
'Q8BGY2', 'O54921', 'P28741', 'P21396', 'P36506', 'P10637',
'Q80ZA5', 'P07895', 'Q60930')
diff_check <- df_1 %>%
filter(Accession %in% acc_2_check)
head(diff_check, 13)
## # A tibble: 13 × 13
## Accession Name `Gene Symbol` `EC log2 fold` `EC est fold` `EC p value`
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 P61751 ARF4_RAT ARF4 -2.21 -4.63 0.0126
## 2 Q64611 CSAD_RAT CSAD -5.30 -39.5 0.0203
## 3 P24268 CATD_RAT CTSD 1.32 2.49 0.0346
## 4 Q641Y8 DDX1_RAT DDX1 -0.622 -1.54 0.0203
## 5 Q641Y8 DDX1_RAT DDX1 -0.159 -1.12 0.758
## 6 Q63622 DLG2_RAT DLG2 1.62 3.08 0.0328
## 7 Q63622 DLG2_RAT DLG2 -0.340 -1.27 0.990
## 8 Q6PDL0 DC1L2_MOUSE DYNC1LI2 1.50 2.82 0.0429
## 9 Q8BGY2 IF5A2_MOUSE EIF5A2 1.10 2.15 0.0446
## 10 Q8BGY2 IF5A2_MOUSE EIF5A2 0.628 1.55 0.0874
## 11 O54921 EXOC2_RAT EXOC2 1.77 3.41 0.0328
## 12 P28741 KIF3A_MOUSE KIF3A 2.23 4.68 0.0429
## 13 P21396 AOFA_RAT MAOA 0.818 1.76 0.0466
## # ℹ 7 more variables: `Coc log2 fold` <dbl>, `Coc est fold` <dbl>,
## # `Coc p value` <dbl>, `interaction log2 effect size` <dbl>,
## # `interaction effect size` <dbl>, `interaction p value` <dbl>,
## # fraction <chr>
tail(diff_check, 13)
## # A tibble: 13 × 13
## Accession Name `Gene Symbol` `EC log2 fold` `EC est fold` `EC p value`
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 P21396 AOFA_RAT MAOA 0.818 1.76 0.0466
## 2 P21396 AOFA_RAT MAOA 0.0829 1.06 0.630
## 3 P21396 AOFA_RAT MAOA 0.127 1.09 0.990
## 4 P36506 MP2K2_RAT MAP2K2 -2.25 -4.74 0.0181
## 5 P10637 TAU_MOUSE MAPT 1.22 2.33 0.0446
## 6 P10637 TAU_MOUSE MAPT 0.120 1.09 0.990
## 7 Q80ZA5 S4A10_RAT SLC4A10 3.02 8.10 0.0211
## 8 P07895 SODM_RAT SOD2 0.914 1.88 0.0487
## 9 P07895 SODM_RAT SOD2 0.851 1.80 0.114
## 10 P07895 SODM_RAT SOD2 0.919 1.89 0.990
## 11 Q60930 VDAC2_MOUSE VDAC2 0.685 1.61 0.0300
## 12 Q60930 VDAC2_MOUSE VDAC2 0.139 1.10 0.386
## 13 Q60930 VDAC2_MOUSE VDAC2 0.0991 1.07 0.990
## # ℹ 7 more variables: `Coc log2 fold` <dbl>, `Coc est fold` <dbl>,
## # `Coc p value` <dbl>, `interaction log2 effect size` <dbl>,
## # `interaction effect size` <dbl>, `interaction p value` <dbl>,
## # fraction <chr>
2025_12_25/26 onwards
Currently only works on SBS Jennifer (insufficient memory on laptop?)
[[fr]][[go]] for nested)Previous Errors
# not needed - works fine with uniprot entries
# # make trimmed entrez_id (rm everything after semi colon)
# df_4 <- df_3 %>%
# ungroup() %>%
# mutate(entrez_id = str_replace(
# GeneID, pattern = ";.*", replacement = ""))
# make vector of fractions (cytolsolic, membrane, & nuclear)
fraction_names <- unique(df_3$fraction)
# GO terms
GO_names <- c('BP', 'CC', 'MF', 'ALL')
# make container list to store results | empty for lists
gsea_GO_results_lists <- vector(mode = "list")
For loop (for each fr fraction),
compute a ranked genelist…
Then nested for loop (for each go
GO term),
if, fr go result is
not empty (> 0) then, make dot plot
else, print message no result.
# works on SBS jennifer
for (fr in fraction_names){ # SEQUENCE - for each item (fr) in fraction_names
# OPERATIONS
# PART 1
# make a ranked gene list for this fraction
geneList_working_df_go <- df_3 %>% # rats only mapped df from Lichti proteomics
# filter by item | fraction is column in df_3, fr is i (current item in seq)
filter(fraction == fr) %>%
# compute ranking metric
mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
# for any duplicate entries, # average & remove any nas
# group_by(Gene.Names..primary.) %>%
group_by(Entry) %>%
summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
# lastly, arrange from largest positive to most negative
arrange(-ranking)
# Part 1.2
# next, deframe into a named vector needed to supply to fgsea()
geneList_vec_working_go <- tibble::deframe(geneList_working_df_go)
geneList_vec_working <- geneList_vec_working_go[!is.na(names(geneList_vec_working_go))]
# store in list
gsea_GO_results_lists$geneList_vec[[fr]] <- geneList_vec_working_go
# Part 1.3 | plots fine but does not save?
# Plot ranked gene List
gL_plot_go <- plot(geneList_vec_working_go, main = paste0(fr, '- Ranked gene list'))
print(gL_plot_go)
gsea_GO_results_lists$geneList_plots[[fr]] <- gL_plot_go
# PART 2
# nested for loop
# so now for each fraction, go through all 3 go terms (BP, CC, MF)
for (go in GO_names){
# operation
# run gseGO for each go term
GO_res_working <- clusterProfiler::gseGO(geneList = geneList_vec_working_go,
OrgDb = org.Rn.eg.db,
keyType = "UNIPROT",
ont = go, #
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 0.05,
pAdjustMethod = "BH"
)
# save GO result original to list vec
gsea_GO_results_lists$GO_results[[fr]][[go]] <- GO_res_working
# extra (not really needed but save anyway
gsea_GO_results_lists$GO_results_tib[[fr]][[go]] <- tibble::as_tibble(GO_res_working)
# PART 3 - if else for dot plot
if (nrow(GO_res_working@result) > 0 ) {
#Make dot plots | this one saves properly
# Part 3.1 - regular dot plots
dp_go <- enrichplot::dotplot(GO_res_working,
showCategory = 30,
font.size = 11,
title = paste0(stringr::str_glue(
'{fr} - {go} - NAc - Lichti 2014 EE vs IC - gseGO'
)))
print(dp_go) # print dot plots
gsea_GO_results_lists$dotplots[[fr]][[go]] <- dp_go # save to list
} else{print(paste0("No GO gsea result for dot plot in ", fr, "-", go, " :("))}
}
}
## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## leading edge analysis...
## done...
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in cytosolic-CC :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in cytosolic-MF :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## leading edge analysis...
## done...
## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 22 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 7 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 1 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 25 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...
## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in membrane-BP :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in membrane-CC :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in membrane-MF :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in membrane-ALL :("
# # horrific name but at least clear...
# Lichti_2014_gsea_GO_results_2025_12_26_v1 <- gsea_GO_results_lists
#
# save(Lichti_2014_gsea_GO_results_2025_12_26_v1,
# file = "output/2025_12_26_Lichti_2014_gsea_GO_results_rats_NAc_v1.rda")
# cow plot for combined + plotlist = to plot a list of plots
# cowplot::plot_grid(plotlist = c(gsea_GO_results_lists$dotplots),
# ncol = 2)
# cow plot not working
# Warning message:
# In as_grob.default(plot) : Cannot convert object of class list into a grob.
# # this works
# print(gsea_GO_results_lists$dotplots)
sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0 LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Pacific/Auckland
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] enrichplot_1.28.4 clusterProfiler_4.16.0 org.Rn.eg.db_3.21.0
## [4] AnnotationDbi_1.70.0 IRanges_2.42.0 S4Vectors_0.46.0
## [7] Biobase_2.68.0 BiocGenerics_0.54.0 generics_0.1.4
## [10] UniProt.ws_2.48.0 gt_1.2.0 lubridate_1.9.4
## [13] forcats_1.0.1 stringr_1.5.2 dplyr_1.1.4
## [16] purrr_1.1.0 readr_2.1.5 tidyr_1.3.1
## [19] tibble_3.3.0 ggplot2_4.0.0 tidyverse_2.0.0
## [22] readxl_1.4.5
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 gson_0.1.0 httr2_1.2.1
## [4] AnVILBase_1.2.0 rlang_1.1.6 magrittr_2.0.4
## [7] DOSE_4.2.0 compiler_4.5.2 RSQLite_2.4.5
## [10] png_0.1-8 vctrs_0.6.5 reshape2_1.4.4
## [13] pkgconfig_2.0.3 crayon_1.5.3 fastmap_1.2.0
## [16] dbplyr_2.5.1 XVector_0.48.0 labeling_0.4.3
## [19] utf8_1.2.6 rmarkdown_2.30 tzdb_0.5.0
## [22] UCSC.utils_1.4.0 bit_4.6.0 xfun_0.53
## [25] cachem_1.1.0 aplot_0.2.9 GenomeInfoDb_1.44.3
## [28] jsonlite_2.0.0 progress_1.2.3 blob_1.2.4
## [31] BiocParallel_1.42.2 parallel_4.5.2 prettyunits_1.2.0
## [34] R6_2.6.1 bslib_0.9.0 stringi_1.8.7
## [37] RColorBrewer_1.1-3 jquerylib_0.1.4 cellranger_1.1.0
## [40] GOSemSim_2.34.0 Rcpp_1.1.0 knitr_1.50
## [43] ggtangle_0.0.9 R.utils_2.13.0 BiocBaseUtils_1.10.0
## [46] igraph_2.1.4 splines_4.5.2 Matrix_1.7-4
## [49] timechange_0.3.0 tidyselect_1.2.1 qvalue_2.40.0
## [52] rstudioapi_0.17.1 yaml_2.3.10 codetools_0.2-20
## [55] curl_7.0.0 rjsoncons_1.3.2 lattice_0.22-7
## [58] plyr_1.8.9 treeio_1.32.0 withr_3.0.2
## [61] KEGGREST_1.48.1 S7_0.2.1 evaluate_1.0.5
## [64] gridGraphics_0.5-1 BiocFileCache_2.16.2 xml2_1.4.0
## [67] Biostrings_2.76.0 ggtree_3.17.0 pillar_1.11.1
## [70] filelock_1.0.3 ggfun_0.2.0 hms_1.1.3
## [73] tidytree_0.4.6 scales_1.4.0 glue_1.8.0
## [76] lazyeval_0.2.2 tools_4.5.2 data.table_1.17.8
## [79] fgsea_1.34.2 fs_1.6.6 fastmatch_1.1-6
## [82] cowplot_1.2.0 grid_4.5.2 ape_5.8-1
## [85] nlme_3.1-168 patchwork_1.3.2 GenomeInfoDbData_1.2.14
## [88] cli_3.6.5 rappdirs_0.3.3 gtable_0.3.6
## [91] R.methodsS3_1.8.2 yulab.utils_0.2.1 sass_0.4.10
## [94] digest_0.6.37 ggrepel_0.9.6 ggplotify_0.1.3
## [97] farver_2.1.2 memoise_2.0.1 htmltools_0.5.8.1
## [100] R.oo_1.27.1 lifecycle_1.0.4 httr_1.4.7
## [103] GO.db_3.21.0 bit64_4.6.0-1